%测试一下三角单元上的积分投影到参考单元上
clc;clear; close all;
f1 = @(x,y) x + y;
f2 = @(x,y) x^2 + y^2;
%int = interior(f1, [2,0],[2,1],[1,0.5]);
%int2 = integral2(f1,1,2,@(x)1-0.5*x,@(x)0.5*x);
node = load("ex_node.mat");
ele = load("ex_ele.mat");
center = load("ex_center.mat");
node = node.ex_node;
ele = ele.ex_ele;
center = center.circucenter;
N = size(ele,1);
M = [1/6,0,0;0,1/6,0;0,0,1/6];
b_value = compute_bvalue();
a_coef = zeros(N,3);
num_value = zeros(N,1);
true_value = zeros(N,1);
%mid_value = compute_mid();
ver_value = compute_vertice();
for i = 1 : N
    tmp_node = [node(ele(i,1)+1,:);node(ele(i,2)+1,:);node(ele(i,3)+1,:)];
    [new_node, tmp_j] = real_and_ref(tmp_node);
    centroid = compute_centroid(tmp_node);
    f_value = zeros(size(new_node,1),1);
    for l = 1 : size(new_node,1)
        f_value(l) = f1(new_node(l,1),new_node(l,2));
    end
    b = zeros(3,1);
    for l = 1 : 3
        fb = f_value.*b_value(l,:)';
        b(l) = compute_int(fb,tmp_j);
    end
    a_coef(i,:) = (tmp_j*M)\b;
    %true_value(i) = f1(node(ele(i,3)+1,1),node(ele(i,3)+1,2));
    true_value(i) = f1(centroid(1),centroid(2));
    %disp(size(a_coef(i,:)));
    num_value(i) = a_coef(i,:) * ver_value;
end
error = sum(abs(true_value-num_value));
%int2 = int2 + integral2(f1,1,2,@(x)x,2);
function int = interior(f, n1, n2, n3)
gauss_points = [0.1666666666667, 0.1666666666667;
  0.6666666666667,0.1666666666667;
   0.1666666666667,0.6666666666667];
 % gauss_points = [1/3 1/3;       % 中心
 %                2/3 1/6;       % 靠近顶点2
 %                1/6 2/3];      % 靠近顶点3
gauss_weights = [0.3333333333333;0.3333333333333;0.3333333333333];
%gauss_weights = [0.1666666666667;0.1666666666667;0.1666666666667];
%计算雅可比矩阵
Jacobi = [-n1(1)+n2(1), n3(1) - n1(1);
    -n1(2)+n2(2), n3(2) - n1(2)];
det_jaco = det(Jacobi);
real_coords = zeros(size(gauss_points)); 
real_vertices = [n1;n2;n3];
for i = 1:size(gauss_points,1)
    xi = gauss_points(i,1);
    eta = gauss_points(i,2);
    
    % 计算形函数
    N1 = 1-xi-eta;
    N2 = xi;
    N3 = eta;
    
    % 计算实际坐标 (x,y)
    x = N1*real_vertices(1,1) + N2*real_vertices(2,1) + N3*real_vertices(3,1);
    y = N1*real_vertices(1,2) + N2*real_vertices(2,2) + N3*real_vertices(3,2);
    real_coords(i,:) = [x, y];
end
integral_gauss = 0;
for i = 1:size(gauss_points,1)
    % 被积函数 f(x,y) = x + y 在实际坐标下的值
    %f = real_coords(i,1) + real_coords(i,2);
    
    % 积分值累加：f * |det(J)| * 权重 * 0.5（三角形面积系数）
    integral_gauss = integral_gauss + f(real_coords(i,1),real_coords(i,2)) * abs(det_jaco) * gauss_weights(i)*0.5;
end
int = integral_gauss;
end

function [new_nodes, det_jaco] = real_and_ref(nodes)
gauss_points = [0.1666666666667, 0.1666666666667;
  0.6666666666667,0.1666666666667;
   0.1666666666667,0.6666666666667];
new_nodes = zeros(size(gauss_points));
Jacobi = [-nodes(1,1)+nodes(2,1), nodes(3,1) - nodes(1,1);
    -nodes(1,2)+nodes(2,2), nodes(3,2) - nodes(1,2)];
det_jaco = abs(det(Jacobi));
for i = 1:size(gauss_points,1)
    xi = gauss_points(i,1);
    eta = gauss_points(i,2);
    
    % 计算形函数
    N1 = 1-xi-eta;
    N2 = xi;
    N3 = eta;
    
    % 计算实际坐标 (x,y)
    x = N1*nodes(1,1) + N2*nodes(2,1) + N3*nodes(3,1);
    y = N1*nodes(1,2) + N2*nodes(2,2) + N3*nodes(3,2);
    new_nodes(i,:) = [x, y];
end
% 求解参数坐标(ξ, η)
%xi = ((x - x1)*(y3 - y1) - (y - y1)*(x3 - x1)) / det_J;
%eta = ((x2 - x1)*(y - y1) - (y2 - y1)*(x - x1)) / det_J;
%mid(1) = ((0.5-nodes(1,1))*(nodes(3,2)-nodes(3,1)) - (0.5-nodes(1,2))*(nodes(3,1)-nodes(1,1))) / det(Jacobi);
%mid(2) = ((nodes(2,1)-nodes(1,1))*(0.5-nodes(3,1)) - (nodes(2,2)-nodes(1,2))*(0.5-nodes(1,1))) / det(Jacobi);
end

function b_value = compute_bvalue()
b1 = @(x,y) 1 - 2.*x;
b2 = @(x,y) 2.*x + 2.*y -1;
b3 = @(x,y) 1 - 2.*y;
gauss_points = [0.1666666666667, 0.1666666666667;
  0.6666666666667,0.1666666666667;
   0.1666666666667,0.6666666666667];
b_value = zeros(3,3);
for i = 1 : 3
    b_value(1,i) = b1(gauss_points(i,1),gauss_points(i,2));
    b_value(2,i) = b2(gauss_points(i,1),gauss_points(i,2));
    b_value(3,i) = b3(gauss_points(i,1),gauss_points(i,2));
end
end

function int = compute_int(f_value,j)
gauss_weights = [0.3333333333333;0.3333333333333;0.3333333333333];
int = f_value' * gauss_weights;
int = 0.5*j*int;
end

function mid_value = compute_mid()
b1 = @(x,y) 1 - 2.*x;
b2 = @(x,y) 2.*x + 2.*y -1;
b3 = @(x,y) 1 - 2.*y;
mid_value = zeros(3,1);
node = [0,0;
    0,1;
    1,0];
centroid = compute_centroid(node);
mid_value(1) = b1(centroid(1),centroid(2));
mid_value(2) = b2(centroid(1),centroid(2));
mid_value(3) = b3(centroid(1),centroid(2));
end

function centroid = compute_centroid(vertices)
    % 输入：vertices为3x2矩阵，每行代表一个顶点的(x,y)坐标
    % 输出：centroid为1x2向量，代表重心坐标(xg, yg)
    
    % 验证输入是否为3个顶点
    if size(vertices, 1) ~= 3 || size(vertices, 2) ~= 2
        error('输入必须是3x2的顶点坐标矩阵');
    end
    
    % 计算重心坐标
    xg = mean(vertices(:, 1));  % x坐标平均值
    yg = mean(vertices(:, 2));  % y坐标平均值
    
    centroid = [xg, yg];
end

function ver_value = compute_vertice()
b1 = @(x,y) 1 - 2.*x;
b2 = @(x,y) 2.*x + 2.*y -1;
b3 = @(x,y) 1 - 2.*y;
ver_value = zeros(3,1);
ver_value(1) = b1(1/3,1/3);
ver_value(2) = b2(1/3,1/3);
ver_value(3) = b3(1/3,1/3);
end